home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / rng / ran1.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-05-14  |  3.0 KB  |  131 lines
  1. /* rng/ran1.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_rng.h>
  23.  
  24. /* This is an implementation of the algorithm used in Numerical
  25.    Recipe's ran1 generator.  It is MINSTD with a 32-element
  26.    shuffle-box. */
  27.  
  28. static inline unsigned long int ran1_get (void *vstate);
  29. static double ran1_get_double (void *vstate);
  30. static void ran1_set (void *state, unsigned long int s);
  31.  
  32. static const long int m = 2147483647, a = 16807, q = 127773, r = 2836;
  33.  
  34. #define N_SHUFFLE 32
  35. #define N_DIV (1 + 2147483646/N_SHUFFLE)
  36.  
  37. typedef struct
  38.   {
  39.     unsigned long int x;
  40.     unsigned long int n;
  41.     unsigned long int shuffle[N_SHUFFLE];
  42.   }
  43. ran1_state_t;
  44.  
  45. static inline unsigned long int
  46. ran1_get (void *vstate)
  47. {
  48.   ran1_state_t *state = (ran1_state_t *) vstate;
  49.  
  50.   const unsigned long int x = state->x;
  51.  
  52.   const long int h = x / q;
  53.   const long int t = a * (x - h * q) - h * r;
  54.  
  55.   if (t < 0)
  56.     {
  57.       state->x = t + m;
  58.     }
  59.   else
  60.     {
  61.       state->x = t;
  62.     }
  63.  
  64.   {
  65.     unsigned long int j = state->n / N_DIV;
  66.     state->n = state->shuffle[j];
  67.     state->shuffle[j] = state->x;
  68.   }
  69.  
  70.   return state->n;
  71. }
  72.  
  73. static double
  74. ran1_get_double (void *vstate)
  75. {
  76.   float x_max = 1 - 1.2e-7f ; /* Numerical Recipes version of 1-FLT_EPS */
  77.  
  78.   float x = ran1_get (vstate) / 2147483647.0f ;
  79.  
  80.   if (x > x_max) 
  81.     return x_max ;
  82.   
  83.   return x ;
  84. }
  85.  
  86.  
  87. static void
  88. ran1_set (void *vstate, unsigned long int s)
  89. {
  90.   ran1_state_t *state = (ran1_state_t *) vstate;
  91.   int i;
  92.  
  93.   if (s == 0)
  94.     s = 1;    /* default seed is 1 */
  95.  
  96.   for (i = 0; i < 8; i++)
  97.     {
  98.       long int h = s / q;
  99.       long int t = a * (s - h * q) - h * r;
  100.       if (t < 0)
  101.     t += m;
  102.       s = t;
  103.     }
  104.  
  105.   for (i = N_SHUFFLE - 1; i >= 0; i--)
  106.     {
  107.       long int h = s / q;
  108.       long int t = a * (s - h * q) - h * r;
  109.       if (t < 0)
  110.     t += m;
  111.       s = t;
  112.       state->shuffle[i] = s;
  113.     }
  114.  
  115.   state->x = s;
  116.   state->n = s;
  117.  
  118.   return;
  119. }
  120.  
  121. static const gsl_rng_type ran1_type =
  122. {"ran1",            /* name */
  123.  2147483646,            /* RAND_MAX */
  124.  1,                /* RAND_MIN */
  125.  sizeof (ran1_state_t),
  126.  &ran1_set,
  127.  &ran1_get,
  128.  &ran1_get_double};
  129.  
  130. const gsl_rng_type *gsl_rng_ran1 = &ran1_type;
  131.